A Periodic Genetic Algorithm with Real- Space Representation for Crystal Structure and 

Polymorph Prediction 
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A novel Genetic Algorithm is described that is suitable for determining the global minimum energy con- 
figurations of crystal structures and which can also be used as a polymorph search technique. This algorithm 
requires no prior assumptions about unit cell size, shape or symmetry, nor about the ionic configuration within 
the unit cell. This therefore enables true ab initio crystal structure and polymorph prediction. Our new algo- 
rithm uses a real-space representation of the population members, and makes use of a novel periodic cut for the 
crossover operation. Results on large Lennard- Jones systems with FCC- and HCP-commensurate cells show 
robust convergence to the bulk structure from a random initial assignment and an ability to successfully dis- 
criminate between competing low enthalpy configurations. Results from an ab initio carbon polymorph search 
show the spontaneous emergence of both Lonsdaleite and graphite like structures. 

PACS numbers: 02.70.-c, 61.50.Ah 



I. INTRODUCTION 

The prediction of crystal structures from first principles has 
long been recognized as one of the outstanding challenges in 
solid state physics 1 ' 2 . The most recent methods of cluster ex- 
pansion assume the lattice structure of the crystal 3,4 . Good 
results for silicon have been shown using minima hopping 5 , 
but this method assumed the number of atoms and unit cell 
of the structures searched. In this communication we demon- 
strate a new method for unbiased ab initio crystal structure 
determination using a novel Genetic Algorithm which makes 
no assumptions of atom number, unit cell or lattice structure. 

Genetic Algorithms (GAs) were first developed by John 
Holland in 19754, and are stochastic global optimization 
methods based on "survival of the fittest". This is a compu- 
tational technique which is used to solve problems in which 
there are many potential solutions, only a small number of 
which are optimal. We initialize our system with a number of 
random candidate solutions which are grouped together into 




FIG. 1: Diagram showing crossover in a fractional representation. 
For each (£, 77) = either (a, b), (b, c) or (c, a) then X = either c, 
a or b. 



a population, with each solution being a member of this pop- 
ulation. There needs to be a way of determining the fitness 
of each member, i.e. some way of telling which members 
are better or worse solutions to the problem than the other 
members. It is this fitness function which defines the prob- 
lem. Population members will be chosen or selected, based 
on their fitness, to become parents to produce offspring in a 
breeding procedure known as crossover. Crossover involves 
creating one or more offspring which are a combination of fea- 
tures from their parents. Each population member needs to be 
encoded or represented in some way such that crossover can 
be performed in a systematic way. For example, in the case 
of binary strings, two parents, 11010011 and 01011001 may 
be split in half and recombined to make two new offspring 
11011001 and 01010011. The offspring may be mutated af- 
ter crossover, which involves making changes to offspring in a 
random way which could introduce beneficial aspects into the 
population. In the binary string case this most often involves 
changing a small percentage of the bits on the string. Using 
each member's fitness the population is updated by only al- 
lowing some population members to survive into the next gen- 
eration, the rest of the population members will be discarded. 
The original formulation used a binary string representation, 
which is described in detail in Holland 6 . In this study the 
problem is the determination of the optimal configuration of 
atoms, where the fitness function is a function of the enthalpy 
of the system. 

GAs were not widely used in the field of solid state physics 
due to the representation issue, the binary nature of which was 
insufficient to describe the complicated atomic and molecu- 
lar systems that are of interest. This changed in 1995 when 
Deaven and Ho 7 described a GA technique for clusters us- 
ing a representation which used the atomic coordinates of the 
systems being studied. Crossover was performed by taking 
a planar cut through the center of each parent and swapping 
halves to generate offspring^. After crossover the resulting 
structures of the resulting population members may be a long 
way from equilibrium and so a quasi-Newton minimizer is 
used to relax each structure to the minimum of the local basin 
of attraction. This reduces the problem by simplifying the po- 
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tential energy surface that is searched by the GA. This method 
shows excellent convergence for a number of complicated 
systems, and recently studies of nanowires 9 and surfacesAQii! 
have been performed using a similar technique. The clus- 
ter calculations are performed using non-periodic codes, but 
the nanowire and surface studies make use of periodic bound- 
ary conditions (PBCs) in either one or two dimensions, and 
a planar cut does not takes these PBCs into account. Bulk 
binary-encoded GA studies have been presented by Wood- 
ley— and real-space encoded studies of bulk molecular crys- 
tals have been reported by Harris et al.^ in which crossover 
was performed by randomly swapping of groups of parame- 
ters between population members, rather than considering the 
structural configuration of the atoms within the crystal. 

In this article, we demonstrate that for periodic GA studies 
using a periodic cut in the crossover operation is superior to a 
planar one. The periodic cut is chosen to have the same peri- 
odicity as the supercell of the population member, and this re- 
duces the discontinuities produced in the offspring operation, 
which would cause extra work for the local minimizer that is 
required in the GA scheme. Results show that a periodic cut 
has a faster convergence than a planar one for large systems. 
We also apply this technique to systems where the unit cell 
of the solution is optimized as well as the atomic coordinates. 
Finally we show that GAs are suitable as a polymorph search 
technique. 



II. METHOD 

The crossover is performed in fractional coordinates, as de- 
scribed in figure H The cut is defined by any periodic func- 
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tion with the same periodicity as the cell, f [s i 

s< atom i s m e fractional position vector for each atom along 
the {(,r)) = (a, b), (b,c) or (c,a) directions. This function 
gives a vector (in fractional coordinates) s cut for each s atom 
in the population member. The metric tensor g = A. T A, where 

h = [a, b, c] (in Cartesian coordinates), is used to calculate 
the product 
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in fractional coordinates (X = c 



when (£, ?/) = (a, b) etc.) and with the criterion 



I > the atom is "above" (outside) the cut 
a cu t \ (2) 
I < the atom is "below" (inside) the cut. 

As the cut is made in fractional coordinates, it does not 
"know" about the Cartesian shape of the cell, so this tech- 
nique allows two population members with different cells to 




FIG. 2: Real-space representation of the periodic cuts in the 
crossover operation. Different wavelengths and amplitudes can be 
used for the cuts along the different cell directions. The cuts are 
calculated in fractional coordinates which allows crossover between 
parents with different cells. The dark gray sections represent one 
part of the cell, the light gray the other, and it is these parts that are 
swapped in crossover. 



be bred during crossover, rather than being constrained to both 
parents having the same supercell. The technique also allows 
the cell size and shape to be evolved along with the crystal 
structure. 

Figure|2]shows the crossover operation where two cuts are 
required in the cell. Every crossover operation performed has 
an equal probability of being calculated with the cuts made in 
reference to either the a, b or c directions. This ensures that 
none of the three co-ordinate directions is preferred over the 
other two, but also allows large areas of each of the population 
members to be undisturbed by the cut. To ensure that no one 
parent is preferred over the other, the center of each of the cuts 
should be made one-quarter and three-quarters up the chosen 
cell axis which gives approximately even mixing. 

In the Deaven and Ho formulation, the plane of the cut was 
defined by the creation of a random unit vector on the sur- 
face of a sphere which was centered on the center of mass of 
the cluster. In our method cuts made in different cell direc- 
tions can have different random wavelengths and amplitudes, 
although a maximum amplitude should be defined so that the 
cut is contained within the simulation supercell, and obviously 
the wavelength of the cut cannot be longer than twice the cell 
vector in that direction, whilst any cut with a wavelength of 
less than half the atomic separation will appear as a flat plane. 

Fitness is determined by the relative enthalpy per atom of 
the population members, and each population member was 
chosen for crossover based on its fitness using roulette wheel 
selection 8 . Similar to Chuang et alJ& the number of atoms in 
each individual population member can be varied by accepting 
all solutions after crossover, or if the number of atoms needs 
to be constrained then solutions are rejected until offspring are 
generated that have the correct number. 

Roulette wheel selection can also be used in the update pro- 
cedure, which determines which members of the original pop- 
ulation should progress through to the next generation. In ad- 
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dition, the lowest energy population member was guaranteed 
to be selected for update when updating by this method. By 
allowing some higher energy members to remain, this update 
procedure prevents the population from becoming stagnated. 
An alternative method for updating is to only allow the low- 
est M population members to proceed to the next generation, 
from a super-population of 2M parents and off-spring. We 
will refer to this method as the "simple" update scheme. 

The way that mutation is performed is also important. This 
is controlled by two quantities, m# and m^. The mutation 
rate is € [0,1], and this determines the probability that 
each atom will be mutated or not after crossover, before the lo- 
cal structure minimization procedure. Once an atom has been 
selected for mutation then it is randomly placed in a cubic 
box with sides of length 2m a which has been centered on the 
atom's original position. 



III. RESULTS 

A. Results from the Empirical Lennard-Jones potential 

All calculations were performed by adding the above GA 
formulation to the ab initio planewave DFT code CASTEP 14 , 
which has also been modified for ease of algorithm testing to 
allow the use of the empirical Lennard-Jones potential^ 
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which has a HCP ground state structure 17 which is almost de- 
generate with the FCC structure (energy difference from HCP 
+0.1 % x k The energy difference between the FCC and HCP 
supercells used in this study was +0.072 %, due to the above 
formulation of the Lennard-Jones potential). The value of a 
was set to 3.405 A, e was set to 120 K, and r cut was set to 
2.5a. While the ground states are very close in energy, to 
switch from FCC to HCP four out of every six layers require 
a stacking fault. 



For the Lennard-Jones results we used a fixed supercell, but 
allowed the number of atoms to either be fixed for each of the 
population members, or be allowed to vary. We compared GA 
minimization calculations using either the planar or periodic 
cuts as described above. The number of population members 
was fixed at M = 16, and the initial number of atoms in each 
population member was set to N = 150 using a hexagonal 
supercell, which is commensurate with perfect FCC and HCP 
structures without stacking faults. The mutation rate, ttir, was 
set to 0.10 (10 %) and the mutation amplitude, itia, was set to 
2. 5 A. The initial configuration of the population members is 
totally randomized, then minimized with the local minimizer 
before proceeding. A total of 200 generations was run for each 
simulation. 



In total 15 simulations were performed from a random start 
for each of the eight combinations of either fixed or variable 
number of atoms, using either the roulette wheel or simple 
update scheme, and with crossover performed using either pe- 
riodic cuts or a planar cut. A summary of the convergence 
times is shown in figure [3] 



1. Empirical Lennard-Jones bulk studies with a fixed number of 
atoms. 



For these studies the number of atoms was kept fixed at 
150 during the whole of the simulation. However, none of the 
60 calculations resulted in minimization down to FCC or HCP 
stacking. The simple update scheme was much faster at reach- 
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FIG. 3: Summary of the convergence times (in number of gener- 
ations) for each variation of the method presented. For each run, 
either periodic or planar cuts were used, using either the roulette or 
the simple update scheme, and the number of atoms could either be 
kept fixed, or be allowed to vary. 
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FIG. 4: Typical results from a 150 atom (variable), 16 population 
member calculation starting from an initial random configuration 
with a mutation rate of 10% and mutation amplitude of 2.5A, us- 
ing roulette wheel selection in the update procedure. Periodic cuts 
were used and the system converged in 29 generations to a structure 
with 150 atoms and a local minimum enthalpy +0.024 % above the 
HCP minimum. 



ing convergence, and using periodic cuts was much faster than 
using a planar cut when using either update scheme. 



2. Empirical Lennard-Jones bulk studies with a variable number 
of atoms. 

Results from Chuang et a/.— showed that allowing the 
number of atoms to vary helped convergence. Our results also 
show this. The use of periodic cuts using the roulette wheel 
update scheme, or the planar cut using the roulette wheel or 
simple update scheme, allowed the system to be minimized to 
a defect-free ground state structure. Periodic cuts were faster 
to convergence than a planar cut, as shown above. 

We found that the system did not converge into a perfect lat- 
tice structure without allowing for variable atom number. Fig- 
ure [4] shows a typical set of results, using periodic cuts with 
roulette wheel selection for update. In this case the system 
converged in 29 generations to the structure shown in figure 
[5] This configuration is an FCC-HCP hybrid with an energy 
difference of +0.024 % from the HCP ground state, which is 
due to a single FCC plane stacking fault. Similar structures 
were also found using a planar cut, but with longer conver- 
gence times. 



tures with four or six atoms per unit cell 20 - 21 . 

Such calculations either require a graph-theoretical enu- 
meration of possible configurations, resulting in an exponen- 
tially growing search space as the number of atoms per unit 
cell increases, or serendipity. By contrast, our GA approach 
is not restricted to any particular form of bonding or chemical 
intuition, and is very efficient at searching high dimensional 
spaces. As an application of our GA to polymorph determi- 
nation we therefore chose to study the four atom per unit cell 
carbon polymorphs. 

All calculations were performed using density functional 
theory (DFT) to treat the electrons for a given configuration 
of ions so no assumptions were made about the nature of the 
underlying bonding. DFT has been shown to be very accurate 
for the calculation of atomic configurations many times - for 
general reviews see Payne et alJ& or Segall et alA For the 
results shown a planewave basis set was used to represent the 

wavefunction, with a 400 eV cutoff and 0.05 A sampling of 
reciprocal space using the Monkhorst-Pack scheme 23 - 24 . Non- 
local ultrasoft psuedopotentials 25 were used to describe the 
electron-ion interaction and the local density approximation 
(LDA) 26 was used to treat exchange correlation effects. The 
LDA was used in preference to the various generalized gradi- 
ent approximation (GGA) functionals, as these are known to 
underbind weakly interacting systems such as graphite sheets. 
The positions of the ions and the unit cell vectors were si- 
multaneously optimized using the quasi-Newton method of 
Pfrommer et alSL, with the reciprocal space sampling den- 
sity and effective planewave cutoff energy maintained at all 
times. 

Only update by roulette wheel selection is used when per- 
forming a polymorph search to allow a large amount of vari- 
ation in the cell bond lengths and angles, and in the ionic po- 
sitions. Out of 8 population members which had all been ran- 
domly initialized, at the end of the tenth generation there were 
five which had a graphite-like character and three which had a 
Lonsdaleite character. A summary of these results, compared 
with diamond, graphite and Lonsdaleite, is shown in table U 
Structure 1 is shown in more detail in figure [5] Within the 
LDA the binding between graphene sheets in graphite is only 
very weakly dependent on the alignment of the sheets, and 
thus any alignment is permissible, as seen in the differences 



Ab initio carbon polymorph studies with a variable 
supercell. 



An interesting alternative application of our GA is in the 
field of polymorph prediction. For example, there has been 
considerable interest recently in carbon polymorphs. System- 
atic ab initio searches been made of sp 3 -hybridized structures 
with four atoms per unit cell 19 and of sp 2 -hybridized struc- 



FIG. 5: Side on view of minimized structure from figure|4] looking 
down the [011] direction. The colors show the mixture of FCC (gray) 
and HCP (black) stacking. 



FIG. 6: Top view of structure 1 of the carbon structures (graphite- 
like) showing the layer stacking. The top layer is colored black and 
the layer below gray for clarity. 

of alignment between figure|S]and figure[7]which shows struc- 
ture 5. Figure [8] shows structure 7 which is Lonsdaleite. Since 
this is a snapshot of the population after 10 generations, any 
configurations of ions are possible, and so diamond structures 
need not be expected. We are currently optimizing the algo- 
rithm and the values of the parameters in order to generate 
as wide a range of structures as possible and results will be 
presented in a forthcoming publication 28 . 



FIG. 8: View of structure 7 of the carbon structures (Lonsdaleite) 
looking down the [110] direction. 

technique, since this does not force the population down into 
a single basin, but allows the search space to be continually 
explored. Future work will include further ab initio studies 
using a variable number of atoms and the application of this 
method to surfaces. A systematic study of the ab initio car- 
bon polymorphs with the genetic algorithm discussed is still in 
progress and will be presented in a forthcoming publication 28 . 



IV. CONCLUSIONS 

We have demonstrated a general method for first principles 
determination of crystal structures using a genetic algorithm 
in a periodic supercell. This technique exploits the inherent 
periodicity in the system in the calculation of crossover be- 
tween parent members of the population by using a periodic 
cut in the crossover operation. This shows much faster conver- 
gence when compared to crossover performed using a planar 
cut, as used in lower dimensional studies. 

There are a number of advantages to this method. First, by 
performing all crossovers in fractional co-ordinates each pop- 
ulation member may be allowed to have unit cells which have 
different sizes and shapes, and if the local minimizer also opti- 
mizes the cell vectors then the optimization process will not be 
biased by choice of initial structure. Indeed, we always start 
from an initial random structure so there is no initial bias to 
any preconceived solution. Secondly, we also suggest that the 
roulette selection method could be used as a polymorph search 




FIG. 7: Top view of structure 5 of the carbon structures (graphite- 
like) showing the mismatch in the layer stacking. The top layer is 
colored black and the layer below gray for clarity. 
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Categorization 


Bond Lengths 


Diamond 




1.52672 A 


1.52672 A 






1.52672 A 


1.52672 A 


Graphite 




1.40731 A 


1.40731 A 


Lonsdaleite 




1.52244 A 


1.52244 A 






1.54583 A 


1.54583 A 


1 


Graphite-like 


1.40696 A 


1.40699 A 


2 


Lonsdaleite-like 


1.52241 A 


1.52241 A 






1.54583 A 


1.54585 A 


3 


Graphite-like 


1.40679 A 


1.40699 A 


4 


Graphite-like 


1.40707 A 


1.40721 A 


5 


Graphite-like 


1.40701 A 


1.40715 A 


6 


Lonsdaleite-like 


1.52225 A 


1.52236 A 






1.54602 A 


1.54629 A 


7 


Lonsdaleite-like 


1.52225 A 


1.52240 A 






1.54569 A 


1.54599 A 


8 


Graphite-like 


1.40736 A 


1.40746 A 



TABLE I: Summary of the eight carbon structures found in poly- 
morph search compared with diamond, graphite and Lonsdaleite. 
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